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1.  INTRODUCTION 


1.1  Project  Overview 

The  work  reported  herein  was  performed  by  the  Georgia  Tech  Research  Institute  (GTRI)  under 
contract  number  W91  INF-13-1-0414,  “Multiple  Time  Series  Node  Synchronization  Utilizing 
Ambient  Reference.”  The  period  of  performance  of  this  effort  goes  from  01  September  2013  to  3 1 
December  2014.  The  GTRI  principal  investigator  was  Dr.  Shean  Phelps,  while  the  sponsor 
technical  contract  was  Dr.  Micheline  Strand  of  the  U.S.  Army  Research  Laboratory’s  Army 
Research  Office  (ARO).  The  program  manager  (PM)  at  the  Defense  Advanced  Research  Projects 
Agency  (DARPA)  was  Dr.  Jim  Gimlett. 


This  grant  was  awarded  under  the  broader  umbrella  of  the  DARPA’ s  Biochronicity  program, 
originally  managed  by  Dr.  Christian  Macedonia.  The  program’s  goal  was  to  advance  the 
fundamental  understanding  of  the  importance  of  timing  and  synchronization  in  complex  biological 
systems.  The  first  phase  of  the  Biochronicity  program  aims  to  identify  the  common  spatio- 
temporal  clock  signatures  in  biological  systems  using  empirically  derived  data.  This  information 
was  then  to  be  used  to  enable  development  of  predictive  algorithms  for  time-dependent  processes. 
The  fundamental  advancements  in  the  understanding  of  timing  in  biological  systems  could  be 
apply  to  many  fields,  one  of  which  is  enhanced  training  of  military  and  non-military  personnel. 


1.2  Project  Summary 

Special  Operators,  like  professional  athletes,  must  maintain  peak  physical  performance  in  and  out 
of  training,  and  like  professional  athletes,  they  could  benefit  from  having  real  time  biometric  data 
used  to  improve  training  and  performance.  In  professional  sports,  the  use  of  wearable  sensor 
networks  may  be  soon  became  commonplace  in  gathering  training  data  for  performance 
evaluation,  but  the  same  is  not  true  for  Special  Operators,  because  of  their  specialized  and 
exclusive  training,  and  the  demanding  operational  conditions.  One  major  requirement  for 
meaningful  data  analysis  and  collective  signal  processing  targeted  to  performance  assessment,  is 
the  need  for  fine  scale  synchronization  among  communicating  nodes  and  across  multiple  domains. 
The  severe  requirements  that  Special  Operators  have  to  obey  do  not  allow  for  common 
synchronization  techniques,  such  as  those  that  rely  on  a  global  clock.  For  this  reason,  techniques 
that  rely  on  a  different  synchronization  mechanism  could  potentially  be  implemented  into  a 
practical  solution  to  provide  real-time  feedback  data  during  training  or  in  field  operations. 

The  research  performed  under  this  seed  grant  approached  the  problem  using  a  combined  hardware 
and  software  approach.  An  analog  sensor  node  with  generic  architecture  was  prototyped  and 
coupled  with  algorithms  designed  to  synchronize  and  classify  multiple  time  series  sampled  at 
possibly  unknown  and/or  variable  rates.  The  sensor  node  works  independently  of  location, 
acquires  multiple  data  modalities,  and  can  be  ruggedized  accordingly  to  the  conditions  of 
operation.  The  algorithms  are  based  on  the  concepts  of  primitives,  for  which  a  complex  activity 
can  be  decomposed  in  a  subset  of  basic  activities.  For  this  work,  a  dictionary  of  basic  activity  was 
derived  using  a  publicly  available  database.  The  problem  of  synchronization  across  different 
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domains  was  approach  from  a  time  perspective  using  resampling  techniques,  and  from  an  event 
perspective  using  correlation-type  measurements. 

During  the  performing  period,  the  following  products  were  developed 


•  A  functional  analog  sensor  node  prototype  with  the  capabilities  to  acquire  acceleration, 
quaternions  and  electrocardiography  signals  (ECG). 

•  A  set  of  tools  for  the  interface  and  the  streaming/download  of  the  sensor  node  data 

•  A  set  of  algorithm  for  the  analysis  of  cross-domain  biometric  data,  consisting  of  the 
following  components: 

-  Preprocessing  and  normalization 

-  Feature  extraction 

-  Statistical  signal  modeling 

-  Validation  and  comparison 


Moreover,  to  facilitate  data  collection  and  hardware  prototyping  the  following  instrumentation  and 
materiel  were  acquired 


•  A  Shimmer  Sensing  laboratory-grade  wearable  data  acquisition  system 

•  A  National  Instruments  high-speed  synchronous  data  acquisition  system 

•  All  parts  and  components  necessary  to  the  realization  of  the  prototype  sensor  node 

The  selection  of  a  high  quality  laboratory-grade  sensor  system  for  the  acquisition  of  multiple 
physiological  parameters  was  important  for  the  characterization  of  the  reference  dictionary  dataset. 
The  research  team  evaluated  several  currently  available  systems  and  based  its  decision  on  several 
factors;  some  of  the  most  important  were  the  following:  1)  Accessibility  to  data;  2)  Flexibility  of 
parameters;  3)  Number  and  type  of  supported  sensors;  4)  Integration  with  other  developing 
environments  (C,  C++,  Matlab,  Java);  5)  Quality,  compactness,  maturity;  6)  Cost. 


During  the  duration  of  the  project,  the  following  datasets  were  acquired  and/or  analyzed 


•  Open  public  PAMPA2  daily  activity  dataset 

•  PhysioNet.org  ECG  activity  dataset  for  normal  conditions 

•  Mixed  activity  dataset  acquired  in  laboratory  settings  using  the  Shimmer  sensors 

All  the  programming,  analysis,  and  algorithm  development  was  done  in  Matlab®,  and  the  sensor 
node  firmware  was  developed  using  MS  Visual  Studio  Tools. 

Because  of  time  constraints  in  the  development  and  implementation  of  the  sensor  node,  no 
controlled  data  useful  for  the  project  was  acquired  using  the  hardware  prototype. 
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2.  MATERIALS  AND  METHODS 


In  this  section,  methods  and  procedures  will  be  described  for  the  algorithm,  hardware  prototyping 
and  dataset  analysis. 


2.1  Approach  and  Introduction 

Human  activity  modeling  and  recognition  using  wearable  sensor  network  is  the  first  step  in 
quantifying  whole  body  motion  during  activity.  New  sensors  have  the  advantage  to  be  smaller, 
energy  efficient  and  more  precise  than  past  generation  sensors.  Moreover,  multiple  sensors  can 
now  be  integrated  in  a  single  package  creating  a  truly  wearable  sensor  suite  that  is  unobtrusive 
enough  to  be  used  during  normal  everyday  activity  or  physical  exercise.  The  nature  of  this  new 
generation  of  devices  makes  them  optimal  tools  for  human  activity  research  studies  in  multiple 
disciplines.  For  instance,  in  healthcare  wearable  sensors  applications  are  found  in  the  study  of 
Parkinson’s  disease,  in  the  monitoring  of  medical  conditions  for  post-traumatic  events,  and  in  the 
study  of  human  motion  for  personalized  treatments.  In  sport  applications,  modem  wearable  sensor 
systems  are  used  to  quantify  athletic  performance  and  for  training  optimization.  In  the  military, 
wearable  systems  are  found  in  multiple  applications,  from  real-time  physiological  heath 
monitoring  to  the  assessment  of  injury  in  the  battlefield. 

The  traditional  approach  in  sensor-based  activity  recognition  considers  activity  as  a  continuous 
motion  and  it  is  based  on  global  features.  Here,  the  data  stream  is  divided  into  windows,  and  each 
window  is  processed  separately.  The  length  of  the  window  is  chosen  to  allow  for  robust  feature 
extraction  and  accurate  classification.  Although  this  model  works  best  in  controlled  setting  and 
laboratory  tests,  real  world  activities  are  characterized  by  high  variability  and  the  collected  signal 
exhibit  a  collective  non-stationarity.  For  this  type  of  signals,  performance  is  tightly  related  to 
window  length,  making  fixed-window  techniques  suboptimal. 

A  different  approach  to  global  features  is  represented  by  the  use  of  statistical  motion  primitive 
models,  where  a  dictionary  is  trained  from  a  set  of  primitive  activities  and  statistical  models  are 
obtained  for  each  set.  This  approach  is  less  sensitive  to  variations  and  can  captures  the  local  aspects 
of  the  activity  signal. 

In  this  work,  the  basic  principle  of  motion  primitive  was  implemented  on  the  statistical  modeling 
of  feature  sets  using  probability  density  functions.  The  primitive  models  were  trained  on  a 
publically  available  machine  learning  dataset  containing  several  representative  samples  of 
everyday  activities.  This  dataset  contained  data  from  sensor  positioned  at  different  location  on  the 
body,  for  18  distinct  activities  (a  complete  description  of  the  dataset  is  presented  later  in  the  text). 
The  data  was  preprocessed  and  representative  features  of  motion  were  extract  using  time, 
frequency  and  statistical  quantities.  Because  of  the  large  size  of  the  feature  matrix,  dimensionality 
reduction  using  principal  component  analysis  was  used.  The  statistical  modeling  was  done  using 
Gaussian  mixture  models  on  the  vector  space  spanned  by  the  principal  components.  A  model  for 
each  activity,  location  and  dataset  available  was  obtained  and  stored  in  a  dedicated  repository. 
After  training,  the  models  were  cross  validated  and  used  for  classification.  Three  use-case 
scenarios  were  developed,  and  they  are  Training,  Evaluation  and  Comparison.  In  Figure  1  is 
illustrated  a  general  diagram  for  the  processing  of  raw  physiological  data,  while  in  Figure  2  is 
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illustrated  the  comparison  of  two  activities.  Explanation  of  the  functionality  of  each  block  is 
presented  later  in  the  text. 


Figure  1. Generic  diagram  illustrating  the  acquisition  and  processing  or  physiological  data.  The  model  derived  from 

the  processing  is  stored  in  a  dedicated  repository. 


> 


-> 


Figure  2.  Generic  diagram  illustrating  the  comparison  of  two  activities  from  two  subjects.  The  result  of  this 
comparison  is  a  chain  of  states  describing  the  activity  time  evolution. 
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2.2  Datasets  used  for  analysis 

Several  datasets  were  used  throughout  the  duration  of  the  project.  Here,  a  brief  summary  for  each 
one  of  them  is  presented. 


PAMAP  2  dataset 

The  UC  Irvine  University  Physical  Activity  Monitoring  for  Aging  People  II  (PAMAP2)  dataset 
was  used  for  the  analysis  and  model  characterization  of  basic  activities  like  sitting,  standing, 
walking,  running  etc.  The  PAMAP  2  dataset  was  used  in  several  peer-reviewed  studies  and 
contains  a  wide  array  of  reference  activities  acquired  for  multiple  subjects.  In  particular,  the 
PAMAP  2  dataset  contains  data  recorded  from  9  subjects  (8  males  and  one  female,  age  27.22  ± 
3.31  years)  and  a  total  of  18  activities.  Each  participant  wore  three  Colibri  wireless  inertial 
measurement  units  (IMU)  sensors  and  one  heart  rate  monitor  (HR-monitor).  Data  was  collected 
at  three  location  on  the  body,  at  the  chest,  at  the  wrist  and  at  the  ankle.  ECG  was  measured  at  the 
chest  only  using  the  HR-monitor.  Sampling  frequency  was  of  100  Hz  for  the  IMUs  and  9  Hz  for 
the  HR-monitor.  The  IMU  collected  acceleration  at  two  different  scales,  +16g  and  ±6 g,  angular 
velocity  ( rad/s ),  magnetometer  (gT)  and  temperature  (°C).  The  HR-monitor  registered  heart  rate 
in  beat-per-minute  (bpm).  For  the  data  collection,  each  participant  followed  a  protocol  containing 
twelve  different  activities.  A  total  of  10  hours  of  data  were  collected,  from  which  nearly  8  hours 
were  labeled  accordingly  to  Table  I. 

The  activities  of  Sitting  (SIT),  Standing  (STD),  Lying  (LNG),  Walking  (WLK),  Running  (RUN)  and 
Cycling  (CYC)  were  chosen  for  the  analysis.  An  example  of  accelerometer  data  from  these 
activities  is  shown  in  Figure  3. 
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Table  1-  PAMPA2  DATASET  ACTIVITIES  DESCRIPTION  (Highlighted  activities  were  selected  for  the  study) 


ID 

Activity 

Description 

ID 

Activity 

Description 

1 

Lying 

Lying  auietlv,  small 

10 

Cycling 

Biking  with  real  bike 

2 

Sitting 

Sitting  in  a  chair 

11 

Running 

Jogging  outside 

3 

Standing 

Standing  still 

12 

Rope  Jumping 

Basic  jumps  and  alternate 

4 

Ironing 

Ironing  1  or  2  T-Shirts 

13 

Watching  TV 

Watch  TV  at  home 

5 

Vacuuming 

Vacuum  cleaning 

14 

Computer  Work 

Normal  office  work 

6 

Ascending  stairs 

Going  upstairs 

15 

Car  Driving 

Driving  normal  traffic 

7 

Descending  stairs 

Going  downstairs 

16 

Folding  Laundry 

Folding  T-Shirts 

8 

Normal  walking 

Walking  at  moderate 

17 

House  Cleaning 

Dusting 

9 

Nordic  walking 

With  walking  poles 

18 

Playing  Soccer 

Running,  dribbling,  passing 

Figure  3.  Example  of  vertical  acceleration  from  chest  mounted  sensor  recorded  for  the  activities  of  Walking, 

Running,  Lying,  Standing,  Sitting  and  Cycling. 


PhysioNet  Datasets 

For  the  analysis  of  physiological  data  from  ECG,  several  PhysioBank1  dataset  were  investigated. 
The  data  from  this  repository  is  often  used  for  peer-reviewed  studies  in  the  biomedical  research 
community  and  it  is  well  documented  and  characterized.  The  datasets  considered  from  this  project 
(listed  below)  were  used  to  derive  the  preprocessing  algorithms  and  the  correct  type  of  features. 


•  AAMI  EC13.  This  dataset  contains  signal  from  10  short  recordings  (60  sec)  specified  by 
the  current  American  National  Standard  for  testing  various  devices  that  measure  heart 


1  www.nhvsionet.erg/nhvsiobank 
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rate.  Each  recording  contains  a  single  lead  ECG  sampled  at  720Hz  and  12  bit  resolution 
over  ±10mV. 

•  MIT-BIH  Arrhythmia.  This  dataset  contains  signal  from  the  first  generally  available  set 
of  standard  test  material  for  the  evaluation  of  arrhythmia  detectors.  The  database  contains 
long  recording  (48  half-hours)  from  two  channel  ambulatory  ECG  from  47  subjects. 
Samples  are  recorded  at  360Hz  with  1 1-bit  resolution  over  ±10mV. 

•  BIDMC  CHF.  This  dataset  contains  long  term  ECG  from  subjects  with  severe  congestive 
heart  failure  (CHF)  conditions.  Individuals  recordings  are  about  20  hours  double  lead  at 
250Hz  at  12-bit  resolution  over  ±10mV.  Original  signal  bandwidth  is  about  0.1  Hz  to  40Hz. 

•  Fantasia.  This  dataset  contains  120-minute  recordings  of  resting  ECG  and  respiratory 
signals  from  48  healthy  subjects,  digitized  at  250Hz. 

An  example  of  ECG  data  from  the  MIT-BIH  and  AAMIecl3  dataset  is  shown  in  Figure  4.  In 
addition  to  ECG,  photoplethysmography  (PPG)  was  also  considered  at  the  beginning  of  the  project, 
but  was  not  included  in  the  final  data  types  for  this  study.  The  PPG  high  sensitivity  to  motion 
artifacts  and  the  difficulties  of  PPG  data  acquisition  outside  of  controlled  settings  made  this 
measurement  not  a  good  candidate  for  motion  classification  and  it  was  the  leading  factor  in  the 
decision  not  to  consider  it  further  into  the  study. 


time  (sec)  time  (sec) 


Figure  4.  Example  of  ECG  signal  from  dataset  MIT-BIH  (left)  and  AAMIecl3  (right). 


Mixed  activity  laboratory  dataset 

This  dataset  was  acquired  using  the  Shimmer  Sensing  wearable  data  acquisition  system.  Because 
the  Shimmer  system  is  modular  and  allows  wireless  communication  between  a  base  station  and 
the  sensors,  a  diversified  set  of  sensors  were  selected  for  the  laboratory  testing.  The  sensor  types 
selected  for  testing  are  accelerometer,  ECG,  and  9DoF  IMU.  Moreover,  this  system  integrates  a 
high  quality  triaxial  accelerometer  in  every  sensor,  eliminating  the  need  for  adding  standalone 
sensors. 
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For  this  dataset,  volunteers  from  GTRI  laboratory  were  instrumented  with  a  chest,  ankle  and  arm 
Shimmer  sensor  and  instructed  to  perform  a  complex  exercise  routine  that  consisted  in  an 
alternated  session  of  walking-running-walking  followed  by  three  jumps.  This  is  illustrated  in  the 
example  of  Figure  5,  where  acceleration  for  the  vertical,  lateral  and  horizontal  direction  is  plotted 
against  time.  Figure  6  instead  shows  an  example  of  ECG  signal  acquired  with  the  Shimmer  sensors 
at  chest  location.  All  data  is  acquired  at  the  sampling  rate  of  512  Hz. 


-30 


0.4  0.6  0.8 


2.2  2.4 


x  10 


Figure  5.  Example  of  vertical,  lateral  and  horizontal  acceleration  recorded  from  a  chest  mounted  Shimmer  unit  for  a 

sequence  of  walk-run- walk  exercise. 


Figure  6.  Example  of  ECG  signal  acquired  with  the  Shimmer  sensors  at  chest  location. 


2.3  Preprocessing 

Signal  preprocessing  represent  the  first  step  in  the  analysis  of  the  input  data  from  the  sensors  and 
it  is  essential  to  the  conditioning  of  the  raw  data  for  the  remaining  processing  steps.  In  the 
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preprocessing  phase,  data  artifacts  are  mitigated  or  eliminated,  and  although  the  applied  techniques 
varied  depending  on  the  signal  type,  some  operations  are  common  to  all  signals.  In  particular,  all 
signals  received  for  analysis  are  processed  to  remove  DC  drifts  and  noise.  Some  specialized 
processing  is  required  for  ECG  signals.  The  following  list  illustrates  the  preprocessing 
implemented  for  each  signal  type: 

•  Acceleration:  detrending  and  denoising 

•  Angular  Velocity:  detrending  and  denoising 

•  ECG:  detrending,  denoising,  and  removal  of  respiration  signal 


Offset  removal  and  detrending 

When  the  trend  is  a  constant,  the  DC  offset  of  the  signal  can  be  corrected  by  subtracting  the  mean 
of  the  signal  from  the  signal  itself.  This  is  done  by  selecting  a  stationary  portion  of  the  signal, 
usually  starting  from  the  first  or  the  last  sample,  and  computing  the  mean  over  a  window.  This  is 
a  good  estimate  of  the  DC  offset,  and  it  is  then  removed  from  the  signal. 

When  signals  are  characterized  by  variable  offsets,  trends  are  estimated  and  removed  using 
interpolation  techniques.  In  this  case,  a  low  order  polynomial  fit  is  computed  on  the  signal  samples 
to  find  the  best  approximation  to  the  data.  Keeping  the  degree  of  the  fitting  polynomial  low  ensures 
that  the  fit  approximates  the  general  trend  of  the  data  without  removing  important  information. 


Denoising 

Because  sensor  signals  will  invariantly  contain  noise,  the  process  of  removing  or  attenuating  this 
noise  is  important  to  assure  that  the  best  quality  possible.  Moreover,  attenuation  of  unwanted  noise 
in  the  data  is  particularly  important  for  the  analysis  of  signals  tightly  correlated  with  the 
phenomenon  under  observation.  This  is  because  of  salient  features  in  the  data  that  may  be  masked 
by  the  noise,  and  preserving  these  features  guarantees  good  overall  modeling  performances. 

Classical  denoising  techniques  based  on  lowpass  filtering  rely  on  the  assumption  that  the  signal 
and  the  noise  are  separated  in  the  frequency  spectrum.  In  particular,  it  is  commonly  assumed  that 
the  signal  is  contained  in  the  lower  portion  of  the  frequency  spectrum  while  the  noise  is  associated 
with  higher  frequency  components.  Although  this  assumption  holds  especially  well  in  the  case  of 
localized  noise  components,  such  as  frequency  peaks  associated  with  resonant  frequencies,  real 
noise  is  usually  found  over  the  entire  Nyquist  interval.  In  particular,  any  disturbance  that  can  be 
model  as  a  sequence  of  uncorrelated  random  variables  with  given  probability  distribution  function 
follows  in  this  category.  For  a  Gaussian  distributed  discrete  random  variable  n(k),  we  have  white 
Gaussian  noise,  which  power  spectrum  is  given  by  the  Fourier  transform  of  its  autocovariance 
function  Cn  (k)  =  cr^d(k),  or  Pn  (/)  =  o}x  V/  £  (co,  co).  Therefore,  if  signal  and  noise  are 
intermixed  in  the  frequency  range  of  interest,  lowpass  filtering  the  signal  is  not  a  good  solution  for 
the  mitigation  of  the  noise. 

An  alternative  approach  is  offered  by  the  application  of  wavelet  theory  to  signal  analysis.  While 
lowpass  filtering  is  a  form  of  denoising  in  the  frequency  domain  alone,  wavelet  denoising  operates 
in  a  different  transformed  domain.  A  wavelet  is  a  basis  function  that  is  concentrated  in  both  time 
and  frequency,  which,  and  unlike  the  Fourier  transform,  provides  a  scale  dependent  time-frequency 
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representation  of  a  signal.  The  discrete  wavelet  transform  (DWT)  algorithm  represents  the  signal 
using  wavelets  as  basis  functions  through  dilations  and  translations  on  a  dyadic  grid. 

Let  us  define  a  discretized  signal  y(/t)  and  suppose  that  there  are  N  samples  of  a  function  /(t), 
such  that 


Vk  =  /(tfc)  +  k  =  1>2, ...  N  (1) 

where  ek  are  lid1  samples  for  a  Gaussian  process  with  zero  mean  and  unit  variance,  and  on 
represent  the  power  of  the  noise.  The  DWT,  yields  the  following 

Yjk  =  wjk  +  <Tnejk  (2) 

where  wJk  are  the  wavelet  coefficient  of  f(tk ). 

The  advantage  of  the  DWT  representation  is  that,  due  to  linearity  of  the  wavelet  transform,  the 
coefficients  of  the  observed  signal  can  themselves  be  considered  a  noisy  version  of  the  wavelet 
coefficients  of  the  original  signal.  Also,  because  of  the  sparsity  properties  of  the  wavelet  transform, 
the  energy  of  the  original  signal  is  usually  concentrated  in  few  coefficients,  with  the  rest  of  them 
being  very  small  or  close  to  zero.  Therefore,  for  a  signal  contaminated  by  Gaussian  noise  the  DWT 
produces  a  small  number  of  coefficients  with  high  amplitude  and  a  larger  number  of  small 
magnitude  coefficients.  The  approach  for  which  each  coefficient  is  compared  with  a  threshold  to 
decide  whether  it  constitutes  a  desirable  part  of  the  original  signal  or  not,  is  called  wavelet 
thresholding.  The  thresholding  extracts  the  significant  coefficients  by  setting  to  zero  the 
coefficients  which  their  absolute  value  is  below  a  certain  threshold  level  A.  An  example  of  wavelet 
denoising  for  an  impulsive  signal  is  shown  in  Figure  7. 


Figure  7  -  Example  of  wavelet  denoising. 


2  iid :  independent  identically  distributed 
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ECG  preprocessing 

The  ECG  signal  records  the  electrical  activity  of  the  heart  using  electrodes  placed  directly  on  the 
skin  of  the  participants.  The  signals  recorded  by  ECG  are  the  voltage  variations  of  the  heart  cells 
action  potentials  during  contractions.  The  aim  of  preprocessing  for  ECG  is  to  improve  the  general 
quality  of  the  signals  for  a  more  accurate  analysis  and  feature  extraction,  enhancing  the 
characteristics  of  the  QRS  complex  without  corrupting  the  rest.  Noise  and  other  spurious 
components  can  distort  the  ECG  recordings  to  such  an  extent  that  measurements  of  the  QRS 
complexes  becomes  unreliable  or  not  possible.  The  three  main  types  of  noise  usually  affecting 
ECG  signals  are 

•  Low  frequency  baseline  wonder  caused  by  respiration  and  body  movements 

•  High  frequency  random  noise  cause  by  electrical  interferences  and  muscular  activity 

•  Random  shifts  in  ECG  signals  caused  by  poor  electrode  contact 

In  literature,  a  classical  approach  for  the  preprocessing  of  ECG  recordings  and  the  detection  of 
QRS  complexes  is  the  removal  of  baseline  wonder  and  high  frequency  noise  using  an  adapting 
filtering  approach  and  classical  denoising  techniques.  In  this  work,  we  chose  to  process  the  ECG 
using  wavelet  denoising  for  the  elimination  of  random  noise  and  Variational  Mode  Decomposition 
(VMD)  for  the  mitigation  of  baseline  wonder.  Moreover,  during  the  investigation,  an  adaptive 
filtering  method  based  on  the  Least-Mean  Square  (LMS)  algorithm  and  using  the  jointly  recorded 
acceleration-ECG  signals  were  implement  as  an  alternative  to  VMD  processing.  The  drawback  of 
this  approach  is  that  need  for  acceleration  to  be  recorded  simultaneously  to  the  ECG,  in  addition 
to  the  long  converging  time  of  the  adaptive  processor.  For  this  reason,  the  VMD  technique  is  the 
preferred  method,  although  both  techniques  are  presented  here. 

Variable  Mode  Decomposition  (VMD) 

Here,  a  summary  of  the  VMD  method  as  it  was  applied  to  the  ECG  signal  used  in  this  work  is 
presented.  Because  of  the  space  required  fully  present  the  VMD  theory  and  framework,  the  reader 
is  referred  to  the  bibliography  for  any  additional  information. 

The  VMD  method  is  based  on  the  concept  of  Intrinsic  Mode  Function  (IMF)  and  it  represents  an 
effort  to  extend  and  formalize  the  well-known  Empirical  Mode  Decomposition  (EMD)  proposed 
by  Huang  et  al.  In  the  EMD  method,  the  IMF  are  amplitude  and  frequency  modulated  functions 
of  the  form  mk(t)  =  Ak(t)cos{(pk(t)\  with  instantaneous  frequency  defined  by  fk(t)  =  0fc'(t). 
Therefore,  for  long  enough  intervals,  mk(t)  can  be  considered  to  be  a  pure  harmonic  signal. 
Although  the  original  EMD  algorithm  lacks  a  formal  mathematical  description  and  it  is  sensitive 
to  noise  and  variations  in  sampling,  it  was  used  in  a  large  number  of  applications  in  several 
engineering,  medicine,  and  technology  fields. 

The  VMD  algorithm  proposes  an  intrinsic  and  adaptive  variational  method  which  minimization 
leads  to  the  decomposition  of  the  signal  into  its  principal  modes.  VMD  determines  the  relevant 
bands  adaptively,  and  estimates  the  corresponding  modes  concurrently,  properly  balancing  the 
errors  between  them.  The  goal  of  the  VMD  method  is  to  decompose  an  input  signal  into  a  discrete 
number  of  modes  that  have  specific  sparsity  properties  while  reproducing  the  input.  Here,  it  is 
assume  each  mode  mk  be  mostly  compact  around  a  center  frequency  fk(t),  which  is  to  be 


16 


determined  along  with  the  decomposition.  The  constraint  variational  problem,  as  described  in  the 
original  work  of  Dragomiretskiy  and  Zosso,  can  be  formulated  as  following 

min 

(mk'f  k) 

such  that  Zfc™ tfc(t)  =  fit).  For  the  solution  of  the  minimization  with  respect  to  mk  and  fk,  the 
reader  is  referred  to  the  original  work  cited  in  the  bibliography. 

The  application  to  the  VMD  algorithm  to  the  ECG  signal  produces  a  multi- function  representation, 
where  different  components  are  representative  of  the  different  frequency  characteristics  of  the 
ECG.  Therefore,  low  frequency  variations  associated  to  respiration  appears  in  lower  modes 
(associated  to  low  center  frequencies),  while  high  frequency  noise  from  muscular  activity  and 
other  sources  appears  in  higher  modes.  It  is  also  possible  to  reconstruct  the  original  signal  fit) 
discarding  unwanted  components,  obtaining  in  this  way  a  simplified  version  of/(t).  This  is 
illustrated  in  Figure  8,  where  the  original  (red)  ECG  is  shown  together  with  the  VMD-filtered  ECG 
and  Figure  9. 


Figure  8.  Example  of  ECG  signal  corrupted  by  low  variational  noise  associated  with  respiration  and  movement. 
Here  the  original  signal  (red)  is  shown  along  with  the  VMD-filter  one,  obtained  considering  all  but  the  first  and  last 

components  of  a  six-mode  decomposition. 
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Figure  9.  VMD  modes  obtained  from  the  signal  from  Figure  8. 


Least-Mean  Square  (LMS)  filtering 

Adaptive  filtering  was  used  to  remove  low  frequency  wander  from  respiration  signal  and  motion 
artifacts  associated  with  body  motion.  Acceleration  from  a  co-located  ECG/accelerometer  was 
used  in  an  adaptive  filter  to  estimate  motion  induced  artifacts  of  the  ECG.  The  LMS  algorithm  was 
used  for  the  adaptive  processor.  The  MLS  method  is  an  iterative  approach  for  the  minimization  of 
the  mean  square  error  (MSE)  between  the  primary  and  a  reference  signal.  An  adaptive  filter  based 
on  the  LMS  algorithm  is  schematically  illustrated  in  Figure  10,  where  the  accelerometer  signal  is 
dynamically  subtracted  from  the  ECG  signal  such  that  the  MSE  error  is  minimized.  For  the  details 
on  the  implementation  of  the  LMS  algorithm,  the  reader  is  referred  to  the  bibliography.  Despite 
the  only  indirect  correlation  of  the  acceleration  with  the  motion  artifact  corrupting  the  ECG,  the 
simplicity  and  relative  stability  of  the  LMS  algorithm  produced  good  acceptable  results,  as  shown 
in  Figure  11. 
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Figure  10.  Schematic  illustration  of  adaptive  filter. 


Figure  11.  Example  of  LMS  processing  for  the  removal  of  low  frequency  wonder.  Top:  original  ECG  signal 
corrupted  by  respiration  noise.  Bottom:  clear  ECG  signal  after  LMS  filtering. 


2.4  Segmentation 

In  case  of  data  comprised  of  blocks  of  basic  or  constant  activities,  the  developed  algorithm 
incorporates  an  option  to  segment  the  data  and  process  these  segments  independently.  This  option 
is  important  mostly  when  generating  primitive  models  without  bias.  Each  segment  should  be 
representative  of  a  portion  of  the  signal  where  variance  is  somehow  constant. 

Manual  segmentation  is  based  on  user  input,  using  custom  designed  interfaces  to  guide  the  user  in 
the  selection  of  the  relevant  segments.  While  this  method  allows  for  precise  localization  of  the 
segment  boundaries  and  it  is  most  useful  when  dealing  with  more  complex  data,  its  usage  is  only 
realistic  in  an  ad-hoc  analysis  or  when  the  number  of  traces  to  analyze  is  small.  During  manual 
segmentation,  the  user  is  prompt  to  choose  if  manual  segmentation  is  needed,  and  then  he  is 
presented  with  an  interface  where  flexible  breakpoints  can  be  overlapped  to  the  time  trace, 
selecting  the  intervals  that  will  define  the  data  segments.  This  process  is  shown  in  Figure  12. 

For  large  datasets,  an  automatic  segmentation  routine  was  developed.  Two  methods  for  the 
automatic  segmentation  of  data  were  investigated,  one  based  on  continuous  Changes  in  Variance 
(CV),  and  one  based  on  the  Cumulative  Sum  Method  (CSM).  The  CV  method  looks  for  changing 
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in  the  variance  of  the  signal  associated  with  a  change  in  physical  activity.  Transition  times  are 
computed  at  the  point  of  variation.  CSM  computes  the  statistics  of  the  signal  for  a  running  window 
and  compares  the  log-likelihood  of  the  distribution  of  the  signal  at  the  right  and  left  of  the  window 
middle  point.  Examples  of  the  results  obtain  using  these  two  techniques  are  shown  in  Figure  13 
and  Figure  14. 


Figure  12.  Detail  of  the  custom  interface  for  the  definition  of  manual  segmentation  breakpoints. 


Figure  13.  Example  of  CS  method  results  on  acceleration  measured  at  the  arm  for  multiple  recorded  physical 
activities.  The  transition  points  in  black  are  manually  highlighted  for  comparison  with  dashed  lines. 


20 


150  200  250  300  350  400  450  500  550  600 

time  (sec) 


Figure  14.  Example  of  CSM  method  results  on  acceleration  measured  at  the  arm  for  multiple  recorded  physical 
activities.  The  transition  points  in  black  are  manually  highlighted  for  comparison  with  dashed  lines. 


2.5  Feature  Extraction 

Features  of  interest  were  derived  from  the  isolated  segments  of  data  and  used  for  classification.  A 
dimensionality  reduction  algorithm  based  on  Principal  Component  Analysis  was  developed  and 
used  on  the  derived  feature  sets.  Two  feature  extraction  algorithms  were  developed,  one  for  the 
acceleration  and  one  for  the  ECG  signals,  both  characterized  by  time  and  frequency  features. 


Acceleration  Features 

Acceleration  features  are  computed  from  recorded  acceleration  for  the  x-,  y-,  and  z-axis  of  the 
sensor  reference  system.  Each  feature  is  computed  on  a  moving  window,  so  that  the  resulting 
feature  is  a  vector  of  the  same  length  of  the  signal  in  input.  The  optimal  length  for  the  observability 
window  was  set  to  three  time  the  sampling  frequency,  as  determine  during  experimentation.  This 
value  was  found  to  be  optimal  for  our  application,  using  the  data  provided  by  the  PAMAP2  dataset. 

The  features  implemented  in  the  algorithm  are  the  following 

•  Movement  Intensity  (MI) 

MI[k ]  =  JaxlW  +  ay[/c]2  +  az[/c]2 

•  Mean 

Mar  =^XM/[*] 

i 

•  Variance 

°m/2  =  _  j  ^(M/[t]  —  Umi)2 
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Skewness 


Kurtosis 


SKmi2 


—  nMi)3 

(V  °a r2) 


AT, 


MI 


Median  of  MI  [i] 

10th  and  90th  percentile  of  MI[i]  ( plOMI  and  p90M7) 

Dominant  frequency,  defined  as  the  maximum  of  the  power  spectral  density  of  MI[i] 
Spectral  Entropy,  defined  as  the  entropy  of  the  normalized  power  spectral  density  (PSD), 
and  compute  as 


H  =  -^PSD  log2 


PSD 


In  addition  to  these  basic  features,  two  additional  features  were  considered: 
•  Angle  between  acceleration  vectors  defined  as 

[tan" 


modulus 


-i  (?±\  _  E 
[aj)  4 


2n 


-  n 


where  (i,j)  =  { x,y,z }  and  i  =£  j 
•  Derivative  of  the  angle  as  computed  above 


ECG  Features 

For  heart  activity,  one  of  the  most  effective  parameter  in  predicting  heart  behavior  is  the  heart  rate 
variability  (HRV),  defined  as  the  variation  of  beat-to-beat  intervals  in  heart  rate.  This  quantity  was 
computed  as  the  time  difference  between  the  peaks  of  QRS  complexes  (also  called  RR  intervals), 
and  provides  an  indirect  measurement  of  the  physiological  behavior  of  the  heart.  Several  indicators 
were  derived  from  this  measure,  where  the  most  important  are  detailed  in  Figure  15.  Moreover, 
instantaneous  heart  rate  (HR)  is  directly  derived  from  the  HRV  signal.  In  literature,  HRV  was  used 
to  quantify  improvement  in  training  and  fitness  performance  in  different  conditions  and  sport 
activity,  and  it  therefore  represent  a  valuable  source  of  information.  Figure  16  to  Figure  18  show 
an  example  from  a  47-minute  long  ECG  recording  using  the  Shimmer  system,  the  derived  HRV, 
HR  and  power  spectrum. 
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Measure 

Description 

Unit 

AVNN 

Average  of  all  NN  intervals 

ms 

SDNN 

Standard  deviation  of  all  NN  intervals. 

ms 

RMSSD 

The  square  root  of  the  mean  of  the  sum  of  the  squares  of 
differences  between  adjacent  NN  intervals 

ms 

pNN50 

Percentage  of  differences  between  adjacent  NN  intervals 
that  are  >  50  ms 

% 

TOTPWR  Total  spectral  power  of  all  NN  intervals  0-0.4  Hz. 

ms2 

VLF 

Total  spectral  power  of  all  NN  intervals  0-0.04  Hz 

ms2 

LF 

Total  spectral  power  of  all  NN  intervals  0.04-0.15  Hz 

ms2 

HF 

LF/HF 

Total  spectral  power  of  all  NN  intervals  0.15-0.4  Hz 

Ratio  of  low  to  high  frequency  power 

ms2 

Figure  15  -  Definition  for  HRV  time  and  frequency  measurements. 
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Figure  16-  (T op)  ECG  signal  from  a  47 -minute  long  recording  using  the  shimmer  system.  (Bottom)  Derived  HRV  measure. 
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Figure  18  -  Power  Spectral  Density  from  the  HRV  signal  of  Figure  16. 


TABLE  2  -  HRV  STATISTICS  COMPUTED  ON  THE  SIGNAL  FROM  FIGURE  16. 


Time  domain  measures 

Frequency  domain  measures 

Mean  HR 

95.4  bpm 

VLF  peak 

9.7  lO"5  Hz 

Std  HR 

12.0  bpm 

VLF  abs  pwr 

0.0029  s2 

AVNN 

0.639  ms 

VLF  rel  pwr 

65% 

SDNN 

0.0784  ms 

LF  peak 

0.04  Hz 

RMSSD 

0.038  ms 

LF  abs  pwr 

0.0011  s2 

LF  rel  pwr 

25% 

HF  peak 

0.15  Hz 

HF  abs  pwr 

4.6 10-4  s2 

HF  rel  pwr 

10% 

LF/HF 

2.38 
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2.6  Statistical  Modeling 

Acceleration  and  ECG  features  distributions  are  modelled  after  dimensionality  reduction,  keeping 
feature  sets  separated  for  identification.  Dimensionality  reduction  is  done  using  Principal 
Component  Analysis  (PCA),  while  the  modeling  of  the  feature  densities  is  done  using  Gaussian 
mixture  models. 

Dimensionality  Reduction  using  PCA 

Principal  Component  Analysis  (PCA)  is  used  to  concentrate  the  information  provided  by  the 
feature  vector  in  few  components  of  high  variance.  The  basic  idea  in  the  application  of  PCA  to  a 
is  to  find  the  component  vectors  ylt  y2, . . . ,  yw  that  explain  the  maximum  amount  of  variance 
possible  by  N  linearly  transformed  components.  An  intuitive  way  to  define  PCA  is  to  use  a 
recursive  formulation.  If  we  define  the  first  principal  component  as  v1,  its  direction  is  found  by 
the  solution  of  the  maximization  problem  vx  =  argmax\\v\\=1  E[(ylX)2],  where  v1  is  the  first 
principal  component  and  X  is  the  vector  containing  the  data.  Therefore,  v1  is  the  projection  of  X 
in  the  direction  that  maximizes  the  variance.  Each  remaining  N  —  1  components  are  computed  by 
repeating  this  process  in  the  remaining  orthogonal  subspace.  The  principal  components  are  then 
given  by 

yt  =  vjx  Vi  =  1,  ...N 

This  operation  transforms  the  column  of  X  into  yt  using  the  quantities  vj .  These  quantities  can 
also  be  computed  using  the  simple  covariance  matrix  C  =  XTX,  where  the  vt  are  the  eigenvectors 
of  C,  which  correspond  to  an  equal  number  of  eigenvalues.  The  transformation  that  from  X  obtains 
yi  using  the  eigenvectors  of  the  covariance  matrix  is  also  known  as  Karhunen-Loeve  transform. 

In  the  analysis,  the  dimensionality  of  the  data  is  reduced  from  N  to  p,  where  p  <  N  to  remove 
unwanted  components.  This  comes  from  the  assumption  that  the  data  in  the  last  N  —  p  components 
is  mostly  noise,  therefore  keeping  the  first  p  components  we  ensure  that  all  signal  information  is 
represented.  For  our  application,  p  is  selected  such  as  to  keep  at  least  80%  of  the  original  signal 
energy.  Although  the  number  of  components  for  which  this  condition  is  met  varies  with  the  type 
of  activity  and  the  measurement  location,  it  was  found  that  on  average  the  first  three  principal 
components  account  for  at  least  80%  of  the  energy.  As  an  example,  Figure  19  shows  the  cumulative 
variance  for  the  activities  of  walking  and  running  measured  at  the  chest  and  ankle. 


Figure  19  -  Cumulative  variance  as  function  of  principal  component,  type  of  activity  and  measurement  location. 
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Gaussian  Mixture  Models 

Each  set  of  principal  components  is  then  modeled  using  Gaussian  mixture  modeling,  a  parametric 
density  estimation  technique.  The  concept  of  mixture  modeling  is  based  on  the  assumption  that 
any  continuous  distribution  can  be  approximated  arbitrarily  well  by  a  finite  mixture  of  normal 
densities  with  common  variance.  For  a  family  of  Gaussian  distributions,  and  a  given  finite  data  set 
X  =  [x1,  x2, ...  xN]  of  N  observations,  the  approximated  mixture  density  has  the  following  form 

K 

/k(x;6»)  =  'YjWkg(x)  9k ) 

fc= 1 

with  0  =  \w1 ...  wK;  9t  02  ...  9K]T  and  9k  =  {fik,  ak}.  The  quantity  wk  is  the  mixing  probability  of 
the  ktb  component  of  the  mixture.  The  model  parameters  wk,  gk  and  ak  need  to  be  estimated  for 
each  one  of  the  K  components  in  the  mixture,  and  the  one-dimensional  Gaussian  kernel  is  defined 
as 

1  G-^fc)2 

gk  =  g(x ;  9k )  =  —  e  2°k2 

\2.7T(7j^ 

A  classical  solution  to  this  problem  is  the  use  of  the  Expectation-Maximization  (EM)  algorithm, 
that  for  a  finite  dataset  and  an  initial  mixture  f(0),  provides  a  sequence  of  mixture  models  f(k^ 
with  increasing  log-likelihood  on  the  data.  For  details  on  the  EM  algorithm,  see  bibliography 
reference. 

The  estimation  of  the  number  of  components  remains  a  problem.  When  a  large  model  order  is 
selected,  the  resulting  model  is  modeling  the  part  of  the  signal  associated  with  noise,  while  on  the 
other  hand  if  a  small  model  order  is  selected,  not  enough  parameters  are  provided  to  the  model  to 
properly  represent  all  signal  features.  In  this  case,  the  model  order  K  was  estimated  as  an  average 
between  repetitions  for  minimum  values  of  the  corrected  Akaike  information  criterion  (AICc)  and 
the  Bayes  Information  Criterion  (BIC).  Here,  the  AICc  is  defined  as 

,  N  2N(3k  - 1) 

AICc  -  -2 £(xn,fk)  +  N_(3k_1)_1 

while  the  BIC  is  defined  as 

BIC  =  -2 l(xn,fk)  +  (3 k  -  1)  log  N 

with  £(xn,fk)  =  £k  =  £m=i  log  fk  as  the  log  likelihood  of  the  estimated  mixture,  and  fk  the  kth 
mixture  component.  An  example  of  computed  mixtures  for  K  =  3  and  two  principal  components 
is  shown  in  Figure  20  for  the  activity  of  walking  and  running. 
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Figure  20  -  Gaussian  mixture  models  of  walking  and  running  activities  estimated  with  two  principal  components. 


2.7  Classification 

The  information  provided  by  the  statistical  model  was  used  to  cluster  data  according  to  its 
characteristic,  the  sensor  location  and  activity  type.  Several  classification  technique  have  been 
investigated  for  this  work,  among  which  Naive  Bayes  classifiers,  and  divergence  metrics 
(Hellinger  distance). 


Naive  Bayes  Classification 

Naive  Bayes  classification  is  an  example  of  supervised  learning  with  the  assumption  that  the 
features  are  independent  for  a  given  class.  Bayesian  classification  is  based  on  the  Bayes  rule  for 
conditional  probability,  and  it  assigns  the  most  likely  class  to  a  give  example.  For  two  events  A 
and  B,  the  Bayes  rule  is  defined  as 


Pr(A\B)  = 


Pr(B\A)Pr(A ) 

Pr(fi) 


A  random  variable  C  is  used  to  denote  the  class  of  an  example,  such  that  C  =  i,  where  ie{  1,2, ...  M } 
and  M  is  the  number  of  models  or  elements  present  in  the  dictionary.  Then,  for  a  given  discriminant 
function  gu  the  classifier  can  be  generally  described  as 
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h(x)  =  arg  ma x  5f(x)  =  {i|V;  *  i:  #y(x)  <  gt{x)} 

For  a  Bayes  classifier,  the  discriminant  function  is  the  class  a-posteriori  probability,  defined  as 


Pi  O )  =  Pr(C  =  i\x)  = 


Pr(x\C  =  i)Pr(C  =  i) 
Pr(x ) 


When  assuming  equally  probable  classes,  Pr(C  =  k)  =  1/M  and  Pr(x)  constant,  the  conditional 
probability  Pr(x|C  =  k)  remains  the  only  term  of  interest.  By  assuming  that  the  features  are 
independent  within  a  class,  the  conditional  probability  can  be  expressed  as  a  product  of  the  single 
conditional  probabilities 

gkNB 0*0  =  Pr(x\C  =  i)Pr(C  =  i)  ~  |~ |Pr(x;|C  =  i)Pr(C  =  i ) 

j 


The  above  expression  implies  that  it  is  possible  to  estimate  the  density  of  each  feature  separately 
and  then  take  the  product  of  the  single  probabilities.  In  our  case,  the  quantity  Pr(x;-  \C  =  i)  is  the 
Gaussian  mixture  associated  with  a  feature  set,  and  using  the  previously  defined  notation 


hNB  (x)  =  arg  max 


Hellinger  Distance 

The  Hellinger  distance  is  divergence  measure  between  probability  distributions.  Assumed  that  P(x) 
and  Q(x)  are  two  probability  distributions  belonging  to  the  same  probability  space,  the  Hellinger 
distance  between  them  is  defined  as 

H(Pj  Q)  =  2  J  dx(/PM-/QM) 

In  this  work,  the  Hellinger  distance  was  used  to  measure  model  divergence  for  the  highlighted 
activities  of  Table  1.  Figure  21  illustrates  the  process  of  estimation  of  the  Hellinger  distance  for 
one-dimensional  and  two-dimensional  Gaussian  mixtures. 


H  ■  D.72D21  -  95fl/D  sigma  bounds  H  =  D.2D59E-  46%  sigma  bounds 


PCI  PCI 

Figure  21  -  Example  of  estimated  Hellinger  distance  for  chest  sensors  CYCLING- WALKING  (left),  and 

RUNNING-SITTING  (right). 
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3.  HARDWARE  PROTOTYPYING 


A  generic  analog  sensor  node  capable  of  acquiring  acceleration  and  ECG  (primary  inputs  of 
interest)  was  developed  as  part  of  this  effort.  The  sensor  is  also  equipped  with  an  IMU  sensor  able 
to  provide  orientation.  The  component  of  the  sensor  node  are  highlighted  in  Figure  22,  and  are  the 
Data  Acquisition,  Processing  and  Transmission  units.  The  Data  Display  and  Analysis  functions 
are  not  included  in  the  sensor  itself  are  part  of  an  external  system.  A  function  block  diagram 
highlighting  the  components  and  the  communication  protocols  used  for  each  functionality  are 
shown  in  Figure  23.  These  are 


•  InvenSense  MPU-9150  -  9  -axis  Motion  Tracking  Inertial  Measurement  Unit  with  3-axis 
MEM  gyro,  3-axis  MEM  accelerometer,  and  3-axis  MEM  magnetometer.  This  unit  has  a 
16-bit  DAC  for  every  sensing  unit. 

•  TI  MSP430  -  16-bit  ultra-low  power  microcontroller,  18MHz  clock. 

•  TI  ADS1293  -  Low  power,  3  channel,  24-bit  analog  front-end  biopotential  measurement 
unit  with  25.6  ksps  maximum  bandwidth 

•  TI  CC2538  -  2.4GHz  IEEE  862.15.4-2006  System-on-a-chip  wireless  microcontroller 


The  firmware  for  the  system  prototype  was  developed  using  the  Contiki3  open  source  OS,  and 
embedded  Real  Time  OS  that  provides  low  overhead,  a  quality  network  stack  and  it  is  optimized 
for  very  low  power  consumption. 


Data  Acquisition/ 
Normalization 


Data  Processing/ 
Transmission 


Data  Display  and 
Analysis 


Figure  22  -  Generic  block  diagram  of  sensor  node  highlighting  major  system  components. 


3  http://www.contiki-os.org/ 
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Figure  23  -  Functional  block  diagram. 


Figure  24  -  Prototype  enclosure  for  the  ECG/IMU  node. 

An  example  of  the  CAD  design  for  the  ECG/IMU  node  is  shown  in  Figure  24.  Based  on  this  CAD, 
a  3D  printed  enclosure  was  manufactured  and  used  to  contain  the  electronics.  Figure  25  shows  the 
CAD  drawing  of  the  circular  ECG/IMU  node  along  with  a  booster  battery  pack  mounted  at 
approximately  center  mass  on  a  shoulder  harness.  Figure  26  shows  a  picture  of  the  sensor  node 
board,  with  standard  ECG  lead  connected  in  black,  red  and  white,  while  Figure  27  shows  an  example 
of  GTRI  researcher  wearing  the  sensor  node  prototype  using  a  shoulder  harness.  This  solution  was 
chosen  as  a  demonstration  only,  and  it  is  not  intended  to  limit  its  usage. 
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Figure  25  -  3D  CAD  drawing  of  the  3D  printed  enclosure  and  a  booster  battery  pack  mounted  at  approximately  at 

the  center  mass  on  a  shoulder  harness. 


Figure  26  -  Picture  detailing  the  sensor  node  main  board.  In  red,  white  and  black  are  visible  the  standard  ECG  lead 

connectors. 
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Figure  27  -  Picture  of  researcher  wearing  the  sensor  node  prototype  and  using  a  commercially  available  harness 

system. 
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4.  RESULTS  AND  DISCUSSION 


4.1  Validation  on  PAMAP2  Dataset 

Data  from  the  PAMAP2  database  was  used  to  build  the  dictionary  of  basic  activity  of  Sitting  (SIT), 
Standing  (STD),  Lying  (LNG),  Walking  (WLK),  Running  (RUN)  and  Cycling  (CYC),  using  the 
algorithm  described  above.  Self- validation  was  performed  on  the  collection  of  model  using  a 
randomly  chosen  datasets  for  each  activity  type  and  for  each  sensor  location.  A  description  of  the 
class  labels  used  in  the  self-validation  is  shown  in  Table  3.  Classification  was  done  using  the 
Hellinger  distance  defined  between  Gaussian  Mixture  models,  and  the  most  likely  label  was 
chosen  as  the  one  that  maximizes  the  Hellinger  distance  over  the  entire  label  set.  Summary  of  the 
average  results  are  shown  in  Figure  28  through  Figure  33.  An  average  classification  rate  is  of  72% 
was  obtained  for  the  cross  validation. 

A  brief  analysis  on  the  results  is  presented  in  the  following: 

•  CYCLYNG  (CYC).  Results  are  shown  in  Figure  28,  where  activity  at  the  wrist  and  ankle 
are  correctly  classified. 

•  LYING  (LYN).  Results  are  shown  in  Figure  29,  where  no  correct  classification  is 
obtained. 

•  RUNNING  (RUN).  Results  are  shown  in  Figure  30,  where  correct  labels  are  assigned  for 
each  location. 

•  SITTING  (SIT).  Results  are  shown  in  Figure  3 1 ,  where  correct  labels  are  assigned  for  the 
chest  and  wrist  locations,  while  the  ankle  location  is  misclassified  as  STD-chest. 

•  STANDING  (STD).  Results  are  shown  in  Figure  32,  where  correct  labels  are  assigned  for 
each  location. 

•  WALKING  (WLK).  Results  are  shown  in  Figure  33Figure  32,  where  correct  labels  are 
assigned  for  each  location. 


TABLE  3-  CLASS  LABELS  USED  FOR  THE  SELF- VALIDATION. 


1 

CYC  ankle 

10 

SIT  ankle 

2 

CYC  chest 

11 

SIT  chest 

3 

CYC  wrist 

12 

SIT  wrist 

4 

LNG  ankle 

13 

STD  ankle 

5 

LNG  chest 

14 

STD  chest 

6 

LNG  wrist 

15 

STD  wrist 

7 

RUN  ankle 

16 

WLK  ankle 

8 

RUN  chest 

17 

WLK  chest 

9 

RUN  wrist 

18 

WLK  wrist 
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Figure  28  -  Self- validation  for  the  activity:  CYCLING.  Correct  labels  are  selected  for  the  activity  recorded  at  the 
wrist  and  ankle,  while  the  chest  data  was  not  correctly  classified. 


Figure  29  -  Self- validation  for  the  activity:  LYING.  No  correct  classification  is  obtained  in  this  case. 


Figure  30  -  Self-validation  for  the  activity:  RUNNING.  Correct  labels  are  assigned  for  each  location. 


Figure  31  -  Self-validation  for  the  activity:  SITTING.  Correct  labels  are  assigned  for  the  chest  and  wrist  locations, 

while  the  ankle  location  is  misclassified  as  STD-chest. 
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Figure  32  -  Self-validation  for  the  activity:  STANDING.  Correct  labels  are  assigned  for  each  location. 


Figure  33  -  Self-validation  for  the  activity:  WALKING.  Correct  labels  are  assigned  for  each  location. 


4.2  Observations  on  laboratory  acquired  data 

Gaussian  modeling  and  classification  was  also  applied  to  data  acquired  in  the  laboratory  for 
controlled  conditions.  For  this  purpose,  the  Shimmer  wearable  sensor  system  was  employed,  and 
data  from  three  locations  on  the  body  were  collected  to  match  the  same  locations  of  the  PAMAP2 
data.  Here,  data  was  acquired  at  a  high  sampling  rate  of  5 12Hz.  The  researchers  wearing  the  system 
were  then  instructed  to  complete  a  routine  that  comprised  of  a  preparation  and  warming  up  session, 
walking  at  a  moderate  pace,  running  at  a  medium  pace  followed  by  a  cooling  down  phase.  The 
exercise  is  always  completed  by  three  jumps,  representing  a  common  denominator  for  the  data. 
This  was  done  because  datasets  were  not  synchronize  among  each  other  and  there  was  the  need  to 
mark  clearly  the  end  of  the  exercise  routine.  Acceleration  and  rotational  velocity  was  acquired  at 
each  node,  and  each  node  acquired  data  independently  of  each  other;  ECG  was  acquired  by  the 
chest  sensor.  An  example  of  acceleration  trace  recorded  at  the  chest  can  be  seen  in  Figure  34,  while 
an  extract  of  low  and  high  intensity  activity  is  shown  in  Figure  35. 

Classification  results  are  shown  in  Figure  36  and  Figure  37.  These  plots  show  the  strength  and 
limitation  of  the  two  approaches.  The  Naive  Bayes  classifier  appears  to  work  better  for  low 
intensity  activities,  which  can  be  compared  to  a  fast  pace  walk  or  slow  run.  In  Figure  36  (left),  the 
highest  score  is  associated  with  class  16,  which  corresponds  to  walking.  The  score  associated  with 
the  sensor  at  the  ankle  does  not  appear  to  contribute  much  to  the  results,  while  the  periodic  motion 
of  chest  and  arm  is  most  likely  driving  the  classifier.  This  is  still  the  case  for  the  higher  intensity 
activity,  which  corresponds  to  a  medium  pace  run,  but  more  uncertainty  can  be  observed,  mostly 
with  the  high  score  of  class  11,12  and  14  corresponding  to  SIT  and  STD. 

Results  for  the  Hellinger  classifier  are  shown  in  Figure  37.  Here,  although  the  selection  of  the  correct 
class  labels  is  not  immediate  due  to  the  high  H  values  across  classes,  decision  on  the  correct  labels 
can  still  be  achieved.  The  maximum  and  second  maximum  of  both  graphs  in  Figure  37  results  in 
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labels  WLK  and  RUN.  This  is  clear  from  the  H  value  associated  with  class  7  and  8,  and  17  for  the 
low  intensity  activity  and  class  8,  17  and  18  for  the  high  intensity  activity. 


Figure  34  -  Example  of  acceleration  trace  recorded  for  a  complex  activity,  recorded  at  the  chest. 
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Figure  35  -  Example  of  acceleration  for  low  intensity  activity  (left)  and  high  intensity  activity  (right) 
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Figure  36  -  Classification  results  for  Naive  Bayes  classifier. 
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Figure  37  -  Classification  results  for  Hellinger  distance  classifier. 


5.  CODE  DEVELOPMENT 

The  algorithm  development  and  code  generation  was  entirely  done  in  Matlab®4,  and  was  organize 
in  three  main  use-case  scenarios.  These  scenarios  are  descried  as  following 


A.  Model  Training:  newly  acquired  activity  data  is  pre-processed,  segmented  and  features 
are  computed.  These  features  are  then  used  to  compute  the  statistical  models  used  later 
for  comparison  and  evaluation  of  new  data. 

B.  Model  Evaluation:  a  newly  acquired  dataset  is  processed  and  its  statistical  model  is  then 
compared  with  stored  models.  A  label  is  assigned  at  the  end  of  this  phase 

C.  Model  Comparison:  two  newly  acquired  dataset  are  processed,  modeled  and  cross- 
compared.  This  procedure  will  provide  a  measurable  alignment  between  two  dataset, 
computed  as  Hellinger  divergence  between  respective  models 


Figure  38  to  Figure  40  show  a  schematic  representation  of  the  use-case  scenarios  described  above.  In 
the  contest  of  the  project,  use-case  A  was  used  to  generate  a  reference  dataset  from  the  PAMPA2 
data,  once  it  is  correctly  preprocessed.  Subsequently,  use-case  scenario  B  was  used  to  evaluate  a 
new  model  towards  the  known  database  of  models,  and  use-case  C  was  used  to  compare  two  model 
between  each  other. 


4  https://www.mathworks.com/ 
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Figure  38  -  Schematic  representation  of  Model  Training  use-case  scenario. 
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Figure  39  -  Schematic  representation  of  Model  Training  use-case  scenario. 


Figure  40  -  -  Schematic  representation  of  Model  Comparison  use-case  scenario. 


During  the  course  of  the  project,  several  interfaces  aimed  at  facilitating  the  interaction  with  the 
code  were  developed.  To  facilitate  the  segmentation  of  complex  datasets,  a  GUI  that  allow  the  user 
to  select  the  desired  breaking  points  was  designed.  The  segmentation  GUI  is  modal,  and  allows 
the  user  to  decide  is  the  segmentation  is  necessary  or  not.  An  example  of  the  modal  GUI  is  shown 
in  Figure  41,  where  the  acceleration  of  the  chest  sensor  is  display  for  allowing  the  user  to  visually 
determine  the  compactness  of  the  data.  If  segmentation  is  chosen,  the  user  proceeds  to  a  second 
GUI,  where  flexible  breakpoints  can  be  overlapped  to  the  time  trace,  selecting  the  intervals  that 
will  define  the  data  segments.  This  process  is  shown  in  Figure  42,  where  a  green  and  red  bar  are 
visible.  When  a  bar  is  first  added,  it  appears  red  in  color,  meaning  that  it  can  still  be  moved  (by 
dragging)  to  the  desired  location.  In  the  moment  the  bar  is  released,  the  software  fixes  its  location, 
recording  it  in  the  list  on  the  right  side  of  Figure  42.  The  segment  then  cannot  be  moved  anymore, 
and  the  user  can  proceed  to  add  a  new  segment  if  desired.  If  a  mistake  is  made,  the  user  can  use 
the  reset  button  and  restart  the  process.  In  this  way,  the  previous  selected  segments  will  be 
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discarded.  An  important  option  is  the  possibility  to  discard  the  first  and  last  segment  of  the  data. 
This  is  often  useful  for  long  dataset,  where  the  data  at  the  beginning  and  the  end  is  often  not 
desired. 


Figure  41  -  Modal  segmentation  GUI. 


Figure  42  -  Segmentation  GUI  details. 


Another  important  addition  to  the  software  is  the  implementation  of  the  Data  Explorer  and 
Labeling  GUI.  This  interface  is  needed  in  the  supervised  learning  process,  where  each  data 
segment  is  labeled  accordingly,  and  the  model  associated  with  it  is  recorded  and  added  to  the 
model  database.  A  view  of  the  Data  Explorer  and  Labeling  GUI  is  shown  in  Figure  43,  where  several 
panels  are  highlighted.  A  description  for  each  panel  is  provided  in  the  following: 
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Panel  1.  This  panel  contains  the  information  about  the  dataset,  and  three  dropdown  menus 
for  the  selection  of  the  trace  of  interest.  By  choosing  the  location  (chest,  arm,  or  ankle),  the 
type  of  signal  (acceleration,  ECG,  and  rotational  velocity)  and  the  segment  number,  the  use 
can  explore  the  entire  dataset.  Only  the  first  channel  of  each  sensor  is  displayed  in  the  plot. 
The  LOAD  button  allows  the  user  to  load  a  new  dataset,  providing  the  option  for  this  GUI  to 
be  used  as  a  stand-alone  product. 

Panel  2.  When  a  Location/Type/Segment  combination  is  selected,  the  software  loads  the 
model  associated  with  the  data,  which  is  previously  computed  and  stored  in  memory.  The 
precise  information  about  the  Gaussian  Mixture  Model  is  displayed  in  Panel  2,  while  the 
model  itself  can  be  display  in  Panel  3. 

Panel  3.  This  panel  displays  the  model  associated  with  the  data  selected  on  Panel  1. 
Models  that  can  be  display  are  one-  and  two-dimensional  models,  while  higher  dimensional 
models  cannot  be  display  in  a  Cartesian  plot.  This  is  because  model  with  more  than  two 
dimension  require  a  four  dimensional  space.  On  the  bottom  of  Panel  3  there  are  the  axis 
zoom  controls  and  a  field  to  control  the  resolution  of  the  model  mesh. 

Panel  4.  This  panel  is  designed  to  allow  the  user  to  easily  label  the  data  selected  on  Panel 
1 .  When  the  user  clicks  on  the  text  box,  the  software  recognize  that,  a  new  label  is  desired 
and  enables  the  element.  Once  a  name  for  the  label  is  chosen,  clicking  the  ASSIGN  Label 
button  creates  an  object  containing  that  model  and  associated  with  that  label.  This  model  is 
save  in  a  dedicated  directory  on  memory.  It  is  also  possible  to  reset  a  label  for  a  given  sensor 
dataset,  by  using  the  RESET  Label  button.  This  process  removes  the  previously  assigned 
label  and  allows  the  used  to  input  a  new  label. 
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Figure  43  -  Data  explorer  and  labeling  GUI  with  highlighted  panels. 


41 


BIBLIOGRAPHY 


P.  Marian  and  T.  A.  Marian,  “Hellinger  distance  as  a  measure  of  Gaussian  discord,”  arXiv: 1408.4477  [ quant-ph /, 
Aug.  2014. 

A.  P.  Dempster,  N.  M.  Laird,  and  D.  B.  Rubin,  “Maximum  Likelihood  from  Incomplete  Data  via  the  EM 
Algorithm,”  Journal  of  the  Royal  Statistical  Society.  Series  B  ( Methodological ),  vol.  39,  no.  1,  pp.  1-38,  Jan.  1977. 

X.  J.  Xie,  J.  J.  Xiang,  and  L.  M.  Liu,  “The  Application  of  Sensor  Technology  in  Physical  Training,”  Applied 
Mechanics  and  Materials ,  vol.  336-338,  pp.  144-147,  Jul.  2013. 

W.  Zijlstra  and  A.  L.  Hof,  “Assessment  of  spatio-temporal  gait  parameters  from  trunk  accelerations  during  human 
walking,”  Gait  &  Posture ,  vol.  18,  no.  2,  pp.  1-10,  Oct.  2003. 

W.  Tao,  T.  Liu,  R.  Zheng,  and  H.  Feng,  “Gait  Analysis  Using  Wearable  Sensors,”  Sensors  (Basel),  vol.  12,  no.  2, 
pp.  2255-2283,  Feb.  2012. 

W.  Zijlstra,  “Assessment  of  spatio-temporal  parameters  during  unconstrained  walking,”  EurJAppl  Physiol ,  vol.  92, 
no.  1-2,  pp.  39-44,  Jun.  2004. 

J.  He,  Y.  Geng,  and  K.  Pahlavan,  “Toward  Accurate  Human  Tracking:  Modeling  Time-of- Arrival  for  Wireless 
Wearable  Sensors  in  Multipath  Environment,”  IEEE  Sensors  Journal ,  vol.  14,  no.  11,  pp.  3996-4006,  Nov.  2014. 

K.  Liu,  T.  Liu,  K.  Shibata,  Y.  Inoue,  and  R.  Zheng,  “Novel  approach  to  ambulatory  assessment  of  human  segmental 
orientation  on  a  wearable  sensor  system,”  Journal  of  Biomechanics ,  vol.  42,  no.  16,  pp.  2747-2752,  Dec.  2009. 

L.  A.  Schwarz,  D.  Mateus,  and  N.  Navab,  “Recognizing  multiple  human  activities  and  tracking  full-body  pose  in 
unconstrained  environments,”  Pattern  Recognition ,  vol.  45,  no.  1,  pp.  1 1-23,  Jan.  2012. 

K.  Altun,  B.  Barshan,  and  O.  Tunsel,  “Comparative  study  on  classifying  human  activities  with  miniature  inertial  and 
magnetic  sensors,”  Pattern  Recognition ,  vol.  43,  no.  10,  pp.  3605-3620,  Oct.  2010. 

J.  J.  Verbeek,  N.  Vlassis,  and  B.  Krose,  “Efficient  Greedy  Learning  of  Gaussian  Mixture  Models,”  Neural 
Computation ,  vol.  15,  no.  2,  pp.  469-485,  Feb.  2003. 

K.  Dragomiretskiy  and  D.  Zosso,  “Variational  Mode  Decomposition,”  IEEE  Transactions  on  Signal  Processing ,  vol. 
62,  no.  3,  pp.  531-544,  Feb.  2014. 

I.  Cleland,  B.  Kikhia,  C.  Nugent,  A.  Boytsov,  J.  Hallberg,  K.  Synnes,  S.  McClean,  and  D.  Finlay,  “Optimal 
Placement  of  Accelerometers  for  the  Detection  of  Everyday  Activities,”  Sensors ,  vol.  13,  no.  7,  pp.  9183-9200,  Jul. 
2013. 

S.  Chernbumroong,  A.  S.  Atkins,  and  H.  Yu,  “Activity  classification  using  a  single  wrist- worn  accelerometer,”  in 
2011  5th  International  Conference  on  Software,  Knowledge  Information,  Industrial  Management  and  Applications 
(SKIMA),  2011,  pp.  1-6. 

L.  Oudre,  J.  Jakubowicz,  P.  Bianchi,  and  C.  Simon,  “Classification  of  Periodic  Activities  Using  the  Wasserstein 
Distance,”  IEEE  Transactions  on  Biomedical  Engineering ,  vol.  59,  no.  6,  pp.  1610-1619,  Jun.  2012. 

J.  Parkka,  M.  Ermes,  P.  Korpipaa,  J.  Mantyjarvi,  J.  Peltola,  and  I.  Korhonen,  “Activity  classification  using  realistic 
data  from  wearable  sensors,”  IEEE  Transactions  on  Information  Technology  in  Biomedicine ,  vol.  10,  no.  1,  pp.  119- 
128,  Jan.  2006. s 

M.  Zhang  and  A.  A.  Sawchuk,  “A  Feature  Selection-based  Framework  for  Human  Activity  Recognition  Using 
Wearable  Multimodal  Sensors,”  in  Proceedings  of  the  6th  International  Conference  on  Body  Area  Networks ,  ICST, 
Brussels,  Belgium,  Belgium,  2011,  pp.  92-98. 


42 


